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Abstrast 

Volatile  organic  compounds  (VOCs)  dissolved  in  the  saturated  zone  are 
transported  into  the  vadose  zone  primarily  by  gaseous  phase  diffusion.  If  the  saturated 
zone  is  remediated,  VOCs  present  in  the  vadose  zone  may  become  a  secondary  source  of 
contamination  for  the  groimdwater.  The  amount  of  VOCs  that  remain  in  the  vadose  zone 
is  dependent  on  site  hydrology,  soil  properties,  and  the  chemical  properties  of  the 
contaminants. 

The  purpose  of  this  study  was  to  determine  what  conditions  caused  VOC 
concentrations  in  the  vadose  zone  to  significantly  recontaminate  the  saturated  zone.  A 
one-dimensional  numerical  model  was  developed  to  investigate  the  transport  of  a  VOC, 
trichloroethylene,  between  the  saturated  and  vadose  zones  under  a  variety  of  conditions. 
The  model  featured  steady-state  imsaturated  water  flow  and  transient  contaminant 
transport.  Transport  mechanisms  included  aqueous  phase  advection-dispersion  and 
gaseous  phase  diffusion.  Partitioning  between  the  water,  gas,  and  soil  compartments 
were  modeled  as  equilibrium  processes.  Sensitivity  analyses  were  performed  on  several 
variables  including  soil  type  (homogeneous  and  heterogeneous  profiles),  water  infiltration 
rate,  and  vadose  zone  depth.  Results  indicated  that  recontamination  was  most  significant 
in  the  presence  of  heterogeneous  soils,  low  infiltration  rates  and  deep  vadose  zones. 
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MASS  TRANSPORT  OF  VOLATILE  ORGANIC  COMPOUNDS 
BETWEEN  THE  SATURATED  AND  VADOSE  ZONES 

L  Introduction 


Background 

Roughly  half  of  the  U.  S.  population  relies  on  groundwater  as  their  primary  source 
of  drinking  water  (Tietenberg,  1994).  As  a  result,  groundwater  contamination  by  organic 
compounds  has  become  a  major  environmental  concern.  Some  primary  sources  of 
contamination  include:  seepage  from  unlined  lagoons  and  other  surface  impoundments; 
leakage  from  pipes,  storage  tanks,  and  equipment;  improper  disposal  practices;  and 
accidental  spills  (Dragun  et  al.,  1984). 

Many  of  these  groundwater  contaminants  are  volatile  organic  compounds  (VOC). 
In  the  saturated  zone,  the  transport  of  dissolved  VOCs  is  normally  dominated  by  aqueous 
phase  advection  and  dispersion.  However,  since  VOCs  have  a  great  tendency  to  partition 
into  the  gaseous  phase,  they  are  readily  transported  into  the  vadose  (or  unsaturated)  zone 
(see  Figure  1).  Once  in  the  vadose  zone,  VOC  transport  is  predominantly  controlled  by 
gaseous  phase  diffusion  (assuming  that  aqueous  phase  advection  is  insignificant  as  a 
transport  mechanism). 
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Figure  1 .  Contamination  of  Vadose  Zone  from  Saturated  Zone 


If  the  groundwater  is  remediated  and  VOC  concentrations  in  the  saturated  zone 
fall  below  those  in  the  vadose  zone,  the  resulting  concentration  gradient  will  cause  VOCs 
to  be  transported  back  into  the  saturated  zone  (see  Figure  2),  Consequently,  VOC  mass 
in  the  vadose  zone  will  become  a  secondary  source  of  contamination  for  the  saturated 
zone.  The  likelihood  that  such  groundwater  recontamination  would  be  significant 
depends  upon  a  myriad  of  factors,  including:  site  hydrology,  soil  properties,  and  the 
physical  and  chemical  properties  of  the  contaminants  themselves. 
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Figure  2.  Recontamination  of  Saturated  Zone  from  Vadose  Zone 


Motivation  for  Research 

Trichloroethylene  (TCE)  and  other  chlorinated  organic  solvents  are  some  of  the 
most  prevalent  organic  contaminants  in  U.  S.  groundwater  supplies  (Dyksen  and  Hess, 
1982;  Conant  et  al.,  1996).  Kerfoot  and  Marrin  (1988)  reported  that  TCE,  a  suspected 
carcinogen,  is  the  most  frequently  identified  substance  at  546  Superfund  sites.  One 
Superfund  site  with  widespread  TCE  contamination  is  McClellan  Air  Force  Base  (AFB), 
California. 
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Nearly  one-third  of  McClellan  AFB  (roughly  four  square  kilometers)  is  underlain 
by  groundwater  that  exceeds  the  maximum  contaminant  level  (MCL)  for  TCE  of  5  j^g/L. 
(Department  of  the  Air  Force,  1994).  Additionally,  the  region  has  experienced  changes  in 
the  direction  of  groundwater  flow  and  a  decline  in  the  water  table  since  the  contamination 
has  occurred.  As  a  result,  the  base  environmental  planners  maintain  that  large  zones  of 
contamination  (principally  TCE)  have  been  created  in  the  subsurface  both  where 
contaminated  groundwater  had  previously  flowed  and  above  where  it  currently  flows.  It 
is  important  to  emphasize  that  this  contamination  is  unlikely  to  contain  nonaqueous  phase 
liquid  (NAPE)  free  product  or  ganglia.  As  far  as  is  known,  contamination  in  these  areas 
was  caused  solely  by  the  dissolved  phase  constituents  in  the  groimdwater.  McClellan’s 
primary  concern  was  that  the  VOCs  that  are  present  in  the  vadose  zone  would 
recontaminate  the  saturated  zone  and  hinder  ongoing  groimdwater  remediation  actions. 
Environmental  planners  speculated  that  restoration  of  these  vadose  zone  areas,  in 
combination  with  saturated  zone  remediation,  might  be  necessary. 

Currently,  the  most  prevalent  vadose  zone  remediation  technology  for  VOC 
contamination  is  soil  vapor  extraction  (SVE)  (Merz  and  Mohr,  1995).  However, 
employing  SVE  to  restore  sites  of  the  magnitude  of  McClellan  AFB  would  be 
prohibitively  expensive.  Additionally,  the  effectiveness  of  SVE  in  these  situations  is 
uncertain  since  the  fate  of  the  contaminants  remains  unknown.  For  example,  if  a 
significant  fraction  of  the  contaminant  mass  either  escaped  into  the  atmosphere  or 
remained  sorbed  to  the  soil  matrix,  the  potential  for  groundwater  recontamination  may  be 
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minimal.  In  such  cases,  remediation  may  be  unnecessary.  Consequently,  in  order  to 
determine  the  need  for  vadose  zone  restoration,  it  is  important  to  understand  the  fate  of 
VOCs  in  the  vadose  zone  and  the  conditions  under  which  significant  groundwater 
recontamination  may  occur. 

Research  Objectives 

(1)  Determine  the  fate  of  VOCs  which  migrate  into  the  vadose  zone  from  the 
saturated  zone. 

(2)  Determine  the  conditions  where  VOCs  in  the  vadose  zone  may  cause 
significant  recontamination  of  groimdwater. 

Thesis  Overview 

The  main  body  of  this  thesis  is  composed  of  five  chapters.  Chapter  I, 
“Introduction,”  provided  a  brief  backgroimd  on  VOC  transport  between  the  saturated  and 
vadose  zones  and  justified  the  focus  of  this  research  in  the  context  of  an  Air  Force  need. 
In  Chapter  II,  “Literature  Review,”  the  works  of  other  authors  are  used  to  identify  the 
predominant  contaminant  transport  mechanisms  and  validate  simplifying  assumptions. 
Chapter  III,  “Methodology,”  describes  the  development  and  validation  of  the  model 
utilized  in  this  investigation.  It  also  provides  a  brief  explanation  of  the  parameters 
selected  for  sensitivity  analysis.  Important  results  are  discussed  in  Chapter  IV,  “Findings 
and  Analysis.”  In  Chapter  V,  “Conclusions,”  the  significant  findings  from  this  study  are 
summarized  and  recommendations  for  future  research  provided.  Additionally,  the 
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attached  appendices  provide  detailed  information  on  notation,  model  output,  equation 
development,  and  model  computer  code. 
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In  the  absence  of  a  NAPL  phase  where  density  effects  and  capillary  forces  may 
influence  contaminant  migration  (Kerfoot  and  Marrin,  1988;  Falta  et  al.,  1989;  Frind  and 


Mendoza,  1990a),  the  potential  transport  mechanisms  controlling  the  movement  of  VOCs 
between  the  saturated  and  vadose  zones  are  limited  to  aqueous  and  gaseous  phase 
diffusion,  dispersion,  and  advection;  and  partitioning  between  the  aqueous,  gaseous,  and 
solid  phases.  The  extent  to  which  each  of  these  processes  contributes  to  the  mass  transfer 
of  VOCs  is  dependent  on  both  the  conditions  in  the  subsurface  and  the  properties  of  the 
contaminant. 

Diffusive  and  Dispersive  Transport 

Diffusion  is  defined  as  the  process  by  which  a  solute  moves  from  an  area  of 
higher  concentration  toward  an  area  of  lower  concentration.  Such  transport  occurs  by 
both  the  random  thermal  motion  of  molecules  (molecular  diffusion)  and  by  their 
interaction  with  pore  walls  (Knudsen  diffusion).  Though  Knudsen  diffusion  has  been 
shown  to  be  significant  in  packed  soil  columns  by  one  study  (Alzaydi  et  al.,  1978),  very 
little  additional  research  on  the  topic  was  attainable.  In  fact,  the  explicit  evaluation  of 
such  diffusion  has  been  generally  ignored  by  soil  scientists  and  hydrologists.  Instead,  its 
effect  is  usually  implicitly  lumped  into  the  tortuosity  factor  (Kreamer  et  al.,  1988). 


Diffilsive  transport  through  the  vadose  zone  occurs  in  both  the  aqueous  and 
gaseous  phases.  However,  since  diffiision  coefficients  for  VOCs  in  air  may  be  several 
orders  of  magnitude  greater  than  those  in  water,  gaseous  diffusion  should  be  the 
predominant  process  (Baehr,  1987;  Sleep  and  Sykes,  1989).  Several  researchers  have 
found  gaseous  diffusion  to  be  a  very  significant  transport  mechanism  (Farmer  et  al., 

1980;  Earp  et  al.,  1982;  Marrin  and  Thompson,  1987).  Nonetheless,  aqueous  phase 
diffusion  may  not  always  be  negligible.  In  unsaturated  soils  possessing  high  water 
contents  (e.g.,  clays)  or  in  areas  near  the  water  table  where  gas  phase  porosity  is  very  low, 
aqueous  diffusion  may  be  an  important  mechanism  for  VOC  transport  (Baehr,  1987). 

Several  authors  have  also  studied  the  effect  of  temperature  on  diffusive  transport 
(Kreamer  et  al,  1988;  Frind  and  Mendoza,  1990b;  Conant  et  al.,  1996).  In  all  these 
studies,  the  researchers  concluded  that  the  vadose  zone  temperature  variations  typically 
encountered  in  most  climates  had  relatively  insignificant  effects  on  diffusive  rates.  Thus, 
it  is  a  reasonable  assumption  to  neglect  the  influence  of  temperature  on  diffusive 
transport. 

Mechanical  dispersion  refers  to  the  mixing  of  a  solute,  as  it  travels  through  a 
porous  medium,  due  to  differential  pore  sizes,  differential  path  lengths,  and  pore  friction 
(Fetter,  1993).  This  mechanism  is  inherently  dependent  on  the  advective  flux  and  it  may 
only  be  neglected  when  advection  is  insignificant. 
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Advection  describes  the  transport  of  a  contaminant  due  to  the  motion  of  the  fluid 
containing  it.  As  with  diffusion,  advective  transport  through  the  vadose  zone  may  occur 


in  both  the  aqueous  and  gaseous  phases.  In  general,  aqueous  advection  is  caused  by 
water  infiltrating  from  the  surface  (McAlary  and  Mendoza,  1990).  Consequently,  this 
transport  mechanism  is  greatly  influenced  by  the  hydraulic  conductivity  of  the  soil  and 
the  precipitation  rate  at  the  surface. 

Gaseous  advection  is  commonly  assumed  to  be  negligible  (Faust,  1985;  Baehr, 
1987;  Abriola,  1989).  This  assumption  is  usually  considered  to  be  reasonable  in  the 
absence  of  a  NAPL  phase  (Abriola  and  Finder,  1986).  Additionally,  Farrier  and 
Massmarm  (1992)  studied  the  influence  that  atmospheric  pressure  variations  induced  on 
advective  transport  in  the  vadose  zone.  They  found  that  though  barometric  pressure  may 
have  notable  short-term  effects,  in  the  long-term,  these  effects  could  be  considered 
negligible.  Also,  others  have  reported  that  water  table  fluctuations  have  had  no 
discernible  effect  on  advective  processes  (Earp  et  al.,  1982;  Kreamer  et  al.,  1988). 

Partitioning 

Partitioning  refers  to  the  process  by  which  a  chemical  becomes  distributed  among 
the  separate  phases  present  in  the  soil  matrix.  Without  a  NAPL  phase,  air/water  and 
soil/water  partitioning  are  usually  the  dominant  mechanisms.  Furthermore,  several 


researchers  have  shown  soil/air  partitioning  (vapor  sorption)  to  only  be  significant  in  very 
dry  soils  (Chiou  and  Shoup,  1985;  Lion  et  al.,  1988;  Culver  et  al.,  1992). 

Additionally,  local  equilibrium  is  commonly  assumed  in  the  literature  (Baehr, 
1987;  Kreamer  et  al,  1988;  Conant  et  al.,  1996).  This  implies  that  “within  some  short 
time  scale  (essentially  instantaneously)  contiguous  phases  reach  a  thermodynamic 
equilibrium.  Thus,  the  mass  fraction  of  a  species  in  one  phase  can  be  related  to  the  mass 
fractions  of  this  same  species  in  the  other  phases  via  partition  expressions”  (Abriola, 
1989).  The  validity  of  this  assumption  at  the  field  scale  has  not  been  confirmed. 
However,  studies  that  employ  the  local  equilibrium  assumption  (LEA)  tend  to 
overestimate  the  extent  of  contamination  (Abriola,  1989),  which  implies  that  it  is  a 
conservative  approach. 

Related  Work 

There  is  a  vast  amoxmt  of  literature  that  concerns  contaminant  transport  in  the 
vadose  zone.  However,  most  of  these  works  focus  on  the  contamination  of  the  saturated 
zone  from  a  NAPL  source  in  the  vadose  zone  (Baehr,  1987;  Falta  et  al.,  1989;  Sleep  and 
Sykes,  1989;  McAlary  and  Mendoza,  1990;  Frind  and  Mendoza,  1990a;  Conant  et  al., 
1996).  Only  a  minority  of  authors  have  studied  the  volatilization  of  contaminants  into  the 
vadose  zone  from  the  saturated  zone 

In  one  field  study.  Barber  et  al.  (1990)  investigated  the  transport  of  methane  from 
the  saturated  zone  to  the  vadose  zone  and  found  that  the  process  was  dominated  by 
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diffusion.  Using  a  simple,  steady-state  diffusion  model,  they  obtained  fair  agreement 
with  field  data,  but  the  assimiption  of  constant  diffiision  with  depth  and  the  neglect  of 
other  processes  such  as  advection  and  sorption  restricts  the  general  applicability  of  the 
model. 

Johnson  and  McCarthy  (1993)  conducted  laboratory  experiments  to  study  the 
transport  of  TCE  from  groundwater  to  the  vadose  zone  and  also  determined  that  diffusion 
was  the  controlling  vertical  transport  mechanism.  They  found  that  TCE  concentrations 
decreased  several  orders  of  magnitude  across  die  capillary  fringe.  Numerical  models 
were  developed  which  incorporated  variable  water  contents  and  diffusion  coefficients 
with  depth.  Both  one-  and  two-dimensional  simulations  agreed  well  with  experimental 
results.  Johnson  and  McCarthy’s  work,  however,  was  limited  to  a  very  shallow 
homogeneous  soil  and  did  not  consider  vertical  water  flow. 

Conclusion 

This  research  advances  our  understanding  of  imsaturated  contaminant  transport. 

It  extends  the  work  of  Johnson  and  McCarthy  (1993)  by  addressing  the  impact  of  soil 
type  (both  homogeneous  and  heterogeneous  soil  profiles),  water  infiltration,  and  vadose 
zone  depth  on  the  transport  of  VOCs  between  the  saturated  and  vadose  zones.  Of 
primary  importance  is  the  potential  for  VOCs  in  the  vadose  zone  to  recontaminate  the 
saturated  zone. 
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III.  Methodology 


Overview 

The  mass  transport  of  the  contaminant  was  estimated  using  a  conservation  of 
mass  approach.  Variable  soil  water  profiles  were  approximated  by  assuming  steady-state 
water  flow.  This  method  yielded  a  partial  differential  equation  which  was  solved  via  a 
finite-difference  approximation.  The  resulting  numeric  model  was  then  validated  against 
a  simplified  analytical  solution.  Sensitivity  analyses  were  conducted  on  soil  type, 
infiltration  rate,  and  vadose  zone  depth. 

Model  Parameters 

For  this  study,  five  different  soil  profiles  were  investigated-three  homogeneous 
and  two  heterogeneous  (see  Figure  3).  In  the  homogeneous  scenarios,  the  vadose  zone 
was  modeled  as  a  single  soil  type  (Sand,  Loam,  or  Clay).  These  soil  types  follow  the 
U.S.  Department  of  Agriculture  (USD  A)  textural  classifications  and  the  parameter  values 
selected  represent  mean  estimates  published  by  Maidment  (1993).  All  model  parameter 
values  are  listed  in  Table  1.  Sand,  loam,  and  clay  were  chosen  to  provide  information 
across  a  wide  range  of  soil  types.  The  heterogeneous  scenarios  included  modeling  a 
sandy  soil  bisected  by  a  50-cm  thick  clay  lens  (Lens),  and  a  system  with  alternating  25- 
cm  layers  of  sand,  loam,  and  clay  (Layered).  The  purpose  of  the  layered  case  was  to 
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Surface 


Sand,  Loam,  or  Clay 


Lens 


Layered 


Figures.  Soil  Type  Profiles 

simulate  very  heterogeneous  conditions.  To  avoid  confusion,  the  five  soil  profiles 
modeled  will  be  categorized  as  simply  soil  types  for  the  remainder  of  this  paper. 

The  aqueous  phase  advective  flux  was  assumed  to  be  constant  and  dominated  by 
the  water  infiltrating  from  the  surface.  The  use  of  a  constant  infiltration  rate  is  most 
reasonable  for  moderate  and  thick  vadose  zones;  it  may  not  be  appropriate  for  porous, 
shallow  soils  whose  water  contents  are  greatly  influenced  by  rainfall  events. 

Nonetheless,  steady  state  infiltration  will  be  assumed  to  be  adequate  for  the  purpose  of 
this  study.  The  infiltration  rates  were  estimated  by  assuming  that  two-thirds  of  the  annual 


precipitation  penetrates  into  the  soil  (T.  Steenhuis,  “Fast  Moving  Solutes  in  the  Vadose 
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Table  1:  Model  Input  Parameters 


USDA  Soil  Classification^ 

Soil  Properties* 

Sand 

Loam 

Clay 

Saturated  volumetric  water  content 

0.417 

0.434 

0.385 

Residual  volumetric  water  content 

0.020 

0.027 

0.090 

Saturated  hydraulic  conductivity  [cm/day] 

500 

30 

1 

Effective  porosity 

0.417 

0.434 

0.385 

Solids  bulk  density  [g/cm^] 

1.55 

1.45 

1.35 

Bubbling  pressure  [cm] 

7.26 

11.15 

37.30 

Pore-size  distribution  index 

0.694 

0.252 

0.165 

TCE  Properties  at  20  °C 

Initial  aqueous  concentration  in  saturated  zone  [mg/L] 

5.0 

Free  water  molecular  diffusion  coefficient  [cm^/day] 

0.729** 

Free  air  molecular  diffusion  coefficient  [cm^/day] 

6993® 

Solids/Water  distribution  coefficient  [cmVg] 

0.118** 

Henry’s  law  constant 

0.350® 

Other  Parameters 

Hydraulic  gradient  (saturated  zone) 

0.02 

Vertical  aqueous  phase  advective  flux  [cm/day] 

0.00 

0.04  0.20 

0.40 

Total  depth  of  vadose  zone  [cm]^ 

300 

1000 

3000 

Longitudinal  dispersivity  in  vadose  zone  [cm]*’® 

30 

100 

300 

Transverse  dispersivity  in  vadose  zone  [cm]** 

3 

10 

30 

“  From  Maidment  (1993) 

Estimated  by  Hayduk  and  Laudie  method  (Lyman  et  a!.,  1982) 


Estimated  by  Fuller,  Giddings,  and  Schettler  method  (Lyman  et  al.,1982) 
^  Estimated  using  regression  from  Chiou  et  al.  (1979) 


From  Howe  et  al.  (1987) 

f 

Values  varied  for  sensitivity  analysis 


®  Values  estimated  using  Peclet  number  of  10 

**  Values  estimated  to  be  1/10  of  longitudinal  dispersivity 


Zone.”  Environmental  Hydrology  Colloquium,  Department  of  Civil  and  Environmental 
Engineering,  University  of  Cincinnati,  OH,  23  February  1996).  Precipitation  rates  were 
selected  to  depict  semi-arid,  temperate,  and  subtropical  climates  (Maidment,  1993).  This 
provided  a  reasonable  range  of  infiltration  rates  for  the  simulations  (see  Table  1).  The 
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zero  infiltration  rate  seenario  served  as  a  baseline  case  from  which  the  reliability  of  the 
model  could  be  verified.  Vadose  zone  thicknesses  vary  widely  in  the  field,  and  the  soil 
depths  modeled  were  chosen  to  reflect  this  attribute.  The  range  of  values  selected  span 
shallow,  moderate,  and  deep  soil  layers  (see  Table  1).  The  VOC  selected  for  this  study 
was  TCE. 

Soil  Water  Profile  Development 

The  water  content  profiles  were  used  to  estimate  the  distribution  of  water  through 
the  entire  vadose  zone  depth.  Unfortunately,  water  content  is  often  discontinuous  in 
space  due  to  varying  matrix  properties  (Abriola,  1989).  This  makes  solving  for  the  water 
content  as  a  function  of  depth,  directly,  very  difficult.  However,  the  pressure  head,  which 
is  always  continuous,  may  be  determined  more  readily. 

Buckingham’s  flux  law,  shown  as  Equation  (1),  describes  the  flow  of  water 
through  an  unsaturated  soil  (Buckingham,  1907): 


(Note:  All  notation  used  in  this  paper  are  defined  in  Appendix  A.) 

The  function  K(T)  in  Equation  (1)  was  expressed  by  the  van  Genuchten  (1980) 
relationship  presented  in  Equation  (2): 

K  Jl  -  (a'E)P-‘  fl  +  (a'F)P 

KQ¥)  =  -L - (2) 

[l  +  (a'E)P]2 
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The  van  Genuchten  soil  parameters  (a,  P,  and  m)  in  Equation  (2)  were  calculated  from 
the  Brooks  and  Corey  parameters  (hb  and  X)  using  the  relationships  shown  in  Equations 
(3),  (4),  and  (5)  (Maidment,  1993).  The  values  used  for  hb  and  X  are  listed  in  Table  1 : 


a  =  (hb)" 

p  =  X  +  l 


m  = 


X.  + 1 


(3) 

(4) 

(5) 


Assuming  that  the  water  infiltration  rate  was  constant,  a  finite-difference 
approximation  of  Equation  (1)  was  solved  iteratively  to  determine  the  pressure  head  as  a 
fimction  of  depth  in  the  vadose  zone,  (4^(z)).  The  solution  technique  began  at  the  water 
table  where,  by  definition,  the  pressure  head  is  zero,  and  propagated  upward  to  the 
surface.  Once  'E(z)  was  determined,  the  water  content  profile  was  calculated  using 
Equation  (6)  (van  Genuchten,  1980): 


0^(z)  =  0,  + 


9s -9, 


(6) 


Soil  water  content  profiles  were  developed  for  each  soil  type,  infiltration  rate,  and 
soil  depth  scenario;  and  are  illustrated  in  Appendix  B.  Also,  Appendix  I  contains  a  more 
detailed  description  of  the  equation  development  used  for  this  model  and  identifies  the 
assumptions  employed.  The  FORTRAN  computer  code  for  the  model  is  listed  in 
Appendix  L. 
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This  model  assumes  that  no  NAPL  phase  exists  and  that  the  air  phase  is  immobile 
(gaseous  advection  negligible).  Additionally,  it  is  assumed  that  the  soil  grains  are 
covered  with  a  continuous  film  of  water;  therefore,  vapor  sorption  is  negligible  (Frind  and 
Mendoza,  1990a).  Since  TCE  does  not  readily  degrade  aerobically,  it  is  further  assumed 
that  all  other  sources  or  sinks  are  insignificant,  as  well.  Consequently,  the  one¬ 
dimensional  mass  transport  equation  reduces  to  Equation  (7); 

where: 


of  water  (upward)  induced  by  the  horizontal  groundwater  flow.  The  parameter,  q^x,  was 
calculated  using  Equation  (1 1)  by  assuming  a  constant  hydraulic  gradient  through  the 
saturated  zone  (the  values  for  6^  and  K(0w)  vary  with  depth  and  were  determined  by  the 


soil  water  content  model).  The  effect  of  this  dispersion  term  only  is  significant  near  the 
water  table  where  the  soil  water  content  (and  hydraulic  conductivity)  is  relatively  high. 


The  other  term,  represents  the  dispersion  due  to  the  infiltrating  water. 


The  contaminant  flux  through  the  porous  media  is  described  by  Equation  (12): 


J  =  (12) 

az  oz 

Assuming  local  equilibrium  and  linear  partitioning,  all  sorbed  and  air 
concentrations  may  be  written  in  terms  of  the  aqueous  concentration  as  shown  in 
Equations  (13)  and  (14)  respectively: 


The  Millington  method  is  used  to  estimate  the  air  and  water  tortuosities  as  shown  in 
Equations  (15)  and  (16)  respectively  (Millington,  1959).  This  method  is  commonly 
employed  in  the  literature  (Baehr,  1987;  Culver  et  al,  1991;  Conant  et  al.,  1996): 


7 


(13) 

(14) 


(15) 


7 


(16) 
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Incorporating  Equations  (13)  through  (16)  into  Equation  (7)  and  assuming  steady  state 
vertical  water  flow,  the  mass  transfer  equation  simplifies  to  Equation  (17): 


^w_ 

<  > 
9wz 

dC^ 

1 

d 

a 

^Rw> 

dz 

Rw 

dz 

where  the  retardation  factor,  R^,  is  described  by  Equation  (18): 

Rw  =  0w  +  +  Pb^-d 

Likewise,  the  contaminant  flux  in  Equation  (12),  reduces  to  Equation  (19): 


J  =  q,,C,-(D,+D,K„)-^ 


(17) 


(18) 


(19) 


Equation  (17)  was  approximated  with  finite  differences  using  a  Crank-Nicolson 
method  modified  to  include  the  advection  term.  The  truncation  error  of  both  of  these 
methods  is  0(Az)^  (order  delta-z  squared).  Equation  (19)  was  also  approximated  as  a 
finite  difference  using  a  backward-difference  method  with  an  error  0(Az)  (Mayers  and 
Morton,  1994).  The  accuracy  of  these  methods  is  assumed  to  be  sufficient  for  the 
purpose  of  this  research.  The  finite-difference  approximations  for  Equations  (17)  and 
(19)  were  solved  using  Dirichlet  boundary  conditions.  The  upper  boundary  (surface)  was 
fixed  at  a  zero  concentration.  To  simulate  the  initial  contamination  of  the  vadose  zone, 
the  lower  boimdary  (water  table)  was  set  at  a  constant  aqueous  concentration  of  5  mg/L, 
which  is  1000  times  greater  than  the  MCL  for  TCE,  and  each  scenario  was  run  for  2000 
simulated  days.  The  resulting  vadose  zone  concentration  profiles  (hereinafter  referred  to 
as  the  initial  concentration  profiles)  are  illustrated  in  Appendix  C.  Then,  to  simulate  the 
recontamination  of  the  saturated  zone,  the  lower  boundary  condition  was  abruptly  set  to 
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zero  and  simulations  continued  until  either  the  highest  aqueous  concentration  in  the 
vadose  zone  fell  below  1  }ig/L  (1/5  of  the  MCL)  or  another  2000  days  was  reached.  A 
more  detailed  equation  development  for  the  contaminant  transport  model  is  included  in 
Appendix  J.  Also,  the  model’s  FORTRAN  computer  code  is  listed  in  Appendix  M, 


Model  Validation 

The  numerical  model  was  validated  against  an  analytical  solution.  If  the  mass 
transport  equation  shown  in  Equation  (17)  is  simplified  by  assuming  that  the  advection 
and  diffusion-dispersion  terms  are  constant  in  space,  the  equation  has  the  analytical 
solution  described  as  Equation  (20); 


where: 


Cw(z,t)  =  -^z  +  e^Dl  J]ak(t)sin 


Uz 


k7iz 


Vk=l 


(20) 


U  = 


0™ 


D  -  +  DgK^j 


(t)=  ^  (l  -  e-"k> ) 

(Ou  ' 


a^(0)  = 


2k7iCi 

_UL  _y_ 

e  cos(k7T) 

f  UL  \ 

e  2®cos(k7i)-l 

V  > 

1 

'  kV^ 

1_  1 

2 

4D^ 

.4D^  j 

(21) 

(22) 

(23) 


(24) 
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(Dk=— +  D^  (26) 

4D 

The  complete  derivation  of  Equations  (20)  through  (26)  is  provided  in  Appendix  K.  A 
comparison  between  the  results  of  the  model  and  the  analytical  solution  shows  that  the 
data  are  nearly  identical  (see  Figure  4).  This  outcome  provides  creditability  for  the  model 
results.  Unfortunately,  no  suitable  field  or  experimental  data  was  attainable  for  further 
verification. 


l.OE-08  l.OE-06  l.OE-04  l.OE-02  l.OE+00 

Relative  Total  Concentration 


Figure  4:  Comparison  of  Model  Results  with  Analytical  Solution 
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Findings  and  Analysis 


Overview 

Sensitivity  analyses  were  conducted  by  varying  soil  type  (both  homogeneous  and 
heterogeneous  profiles),  water  infiltration  rate,  and  vadose  zone  depth.  Mass  balance  was 
checked  by  comparing  the  initial  total  mass  of  contaminant  present  in  the  vadose  zone  to 
the  amount  of  mass  transferred  into  the  atmosphere,  plus  the  amount  remaining  in  the 
vadose  zone,  plus  the  amovint  transferred  back  into  the  saturated  zone.  The  mass  present 
in  the  vadose  zone  was  calculated  by  integrating  the  initial  and  final  concentration 
profiles  over  space.  The  amount  transferred  into  the  atmosphere  or  saturated  zone  was 
determined  by  integrating  the  respective  contaminant  flux  over  time.  Results  which  were 
within  ±10%  of  the  initial  mass  were  considered  to  be  acceptable. 

Unfortunately,  several  scenarios  exhibited  mass  balance  problems  that  could  not 
be  overcome.  Cases  where  high  infiltration  rates  were  coupled  with  deep  vadose  zones 
proved  to  be  most  troublesome.  These  scenarios  experienced  very  rapid  concentration 
drops  which  created  large  errors  in  the  numerical  solution  technique.  The  homogeneous 
clay  and  heterogeneous  lens  soil  profiles  were  also  difficult  to  determine.  Despite  these 
deficiencies,  the  quantity  of  acceptable  results  was  sufficient  to  observe  trends  in  the  data 
and  to  develop  reasonable  conclusions.  Tables  2  and  3  list  the  model  results  for  the 
homogeneous  and  heterogeneous  scenarios  respectively. 
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Table  2:  Model  Results—Homogeneous  Soil  Scenarios 


Sand 


Infiltration  Rate  (cm/day) 

Depth 

0.00 

0.04 

0.20 

0.40 

(m) 

Mi 

1.09 

1.52 

1.80 

1.80 

3 

Mg 

0.70 

0.91 

0.84 

0.61 

Ms 

0.00 

0.00 

0.00 

0.00 

Mw 

0.37 

0.59 

0.93 

1.14 

%Mi 

98 

99 

98 

98 

Mi 

msM 

8.51 

7.72 

5.71 

Mg 

3.24 

3.36 

1.88 

Ms 

0.00 

Mw 

3.14 

5.12 

5.74 

4.98 

%Mi 

99 

97 

Mi 

23.58 

22.37 

7.28 

30 

Mg 

8.30 

3.75 

0.01 

Ms 

0.86 

4.11 

0.25 

Mw 

14.26 

14.34 

7.39 

%Mi 

99 

99 

105 

Clay 

Infiltration  Rate  (cm/day) 

Depth 

0.00 

0.04 

0.20 

0.40 

(m) 

Mi 

0.85 

0.69 

0.57 

E 

3 

Mg 

0.04 

0.00 

0.00 

E 

Ms 

0.34 

0.02 

0.00 

E 

Mw 

0.43 

0.61 

0.54 

E 

%Mi 

95 

91 

95 

E 

All  results  for  clay  with  vadose  zone  depths  of 
10  and  30  m  exhibited  significant  mass  balance 
errors 


Loam 


Infiltration  Rate  (cm/day) 

Depth 

0.00 

0.04 

0.20 

0.40 

(m) 

Mi 

0.77 

1.69 

0.69 

0.59 

3 

Mg 

0.49 

0.47 

0.00 

0.00 

Ms 

0.00 

0.09 

0.00 

0.00 

Mw 

0.26 

1.10 

0.67 

0.58 

%Mi 

98 

98 

97 

99 

Mi 

3.06 

3.49 

1.82 

1.76 

Mg 

1.79 

0.01 

0.00 

0.00 

Ms 

0.35 

0.99 

0.00 

Mw 

0.78 

2.29 

1.66 

1.64 

%Mi 

95 

95 

92 

94 

Mi 

8.63 

5.42 

30 

Mg 

1.04 

0.00 

Ms 

4.43 

1.34 

Mw 

2.71 

3.53 

%Mi 

95 

90 

Notation; 

-  result  exhibited  mass  balance  errors 

Mj  initial  mass  in  vadose  zone 
Mg  mass  transferred  to  the  atmosphere 
Mg  mass  remaining  in  vadose  zone 
M,y  mass  transferred  to  the  saturated  zone 
%  Mj  %  of  Mj  accounted  for  by  mass  balance 
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Table  3:  Model  Results— Heterogeneous  Soil  Scenarios 


Lens _ 

Infiltration  Rate  (cm/day)  Depth 
0.00  0.04  0.20  0.40  (m) 

Mj  3.08  3.00  2.73  -  3 


Ma  0.38  0.09  0.06  — - 

Ms  0.00  0.00  0.00  - 

2.67  2.88  2.73  — - 

■%1^ . '99 . 99 . 102 . = . 

Mj  8.67  10.00  - . -  10 

Ma  1-81  0.44  . - 

Ms  0.00  0.00  - . 

M^  6.83  9.48  . - 

. loi) . 99 . ---- . 

M;  24.94  23.15  . -  30 

Mg  4.39  0.64  . 

Ms  1-15  4.71  . 

M^  18.96  18.66  - . - 

. 98 . 104 . ---- . 


Notation: 

-  result  exhibited  mass  balance  errors 

M;  initial  mass  in  vadose  zone 
Mj,  mass  transferred  to  the  atmosphere 


Layered 


Infiltration  Rate  (cm/day) 

Depth 

0.00 

0.04 

0.20 

0.40 

(m) 

Mi 

2.45 

2.10 

1.27 

1.13 

3 

Ma 

0.65 

0.03 

0.00 

0.00 

Ms 

0.10 

0.10 

0.00 

0.00 

Mw 

1.68 

1.94 

1.23 

1.22 

%Mi 

99 

99 

98 

108 

Mi 

3.23 

3.62 

2.79 

2.65 

10 

Ma 

0.28 

0.00 

0.00 

0.00 

Ms 

0.97 

0.56 

0.06 

0.03 

Mw 

1.92 

2.98 

2.53 

2.49 

%Mi 

98 

98 

93 

95 

Mj 

4.05 

5.62 

30 

Ma 

0.00 

0.00 

Ms 

1.45 

1.09 

Mw 

2.57 

3.99 

%Mi 

99 

90 

— 

Ms  mass  remaining  in  vadose  zone 
M^  mass  transferred  to  the  saturated  zone 
%  Mj  %  of  Mj  accoimted  for  by  mass  balance 


Effect  of  Soil  Type 

Model  results  show  that  soil  type  does  significantly  affect  the  overall  fate  of 
VOCs  in  the  vadose  zone,  particularly  at  deeper  soil  depths.  Moreover,  it  is  apparent  that 
in  certain  instances,  heterogeneities  in  the  soil  profile  may  substantially  increase  the 
potential  for  contaminant  migration  back  into  the  saturated  zone.  Though  each  scenario 
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had  somewhat  unique  results,  Figure  5  illustrates  one  fairly  representative  case.  All  other 
soil  type  results  are  summarized  in  Appendix  D, 


□  Atmosphere 

□  Vadose  Zone 
B  Sat.  Zone 


Soil  Type 


Figure  5:  Effect  of  Soil  Type  on  Contaminant  Fate 
(Infiltration  Rate:  0.04  cm/day;  Vadose  Zone  Depth:  3  m) 


Figure  5  indicates  that,  with  the  exception  of  sand,  only  a  small  fraction  of 
contaminant  mass  is  lost  into  the  atmosphere  and  very  little,  if  any,  remains  in  the  vadose 
zone.  The  majority  of  the  mass  is  transported  back  into  the  saturated  zone.  On  one 
extreme,  thin  sandy  soils  allow  rapid  transport  to  the  surface  and  experience  relatively 
large  losses  to  the  atmosphere.  Conversely,  the  low  conductivity  of  clayey  soil  impedes 
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significant  transport  of  contaminants  from  the  saturated  zone  to  the  vadose  zone.  Thus 
for  both  of  these  cases,  little  contaminant  is  available  for  recontamination  of  the 
groundwater. 

Figure  5  also  shows,  however,  that  soil  heterogeneities  may  create  opportunities 
for  substantial  recontamination  of  the  saturated  zone.  In  these  cases,  highly  permeable 
layers  (e.g.  sand)  encourage  significant  migration  up  through  the  vadose  zone  imtil  a  low 
permeable  layer  (e.g.  clay)  is  encountered  which  hinders  transport.  Thus,  below  a  clay 
layer,  substantial  contaminant  mass  may  accumulate  (see  Figure  6). 


Relative  Total  Concentration 


Figure  6;  Initial  Concentration  Profile  for  a  Clay  Lens  Scenario 
(Soil  Type:  Lens;  Infiltration  Rate:  0.04  cm/day;  Vadose  Zone  Depth:  3  m) 


26 


Figure  6  illustrates  the  initial  contaminant  concentration  profile  for  one  of  the  clay 
lens  scenarios.  It  clearly  shows  a  large  amount  of  mass  present  immediately  below  the 
clay  layer.  Once  the  saturated  zone  is  remediated  and  the  concentration  gradient  is 
reversed,  this  large  mass  will  be  transferred  back  into  the  groimdwater  potentially  causing 
significant  recontamination.  Of  the  two  heterogeneous  soils  modeled,  the  clay  lens 
scenarios  consistently  proved  to  be  the  worst  cases  (returned  most  mass  to  the  saturated 
zone). 

Effect  of  Infiltration  Rate 

Water  infiltration  also  seems  to  have  a  significant  effect  on  the  extent  of 
recontamination.  For  most  soils  imder  normal  infiltration  rates  (between  0.04  and  0.4 
cm/day),  the  mass  returned  to  the  saturated  zone  decreases  with  increasing  infiltration  rate 
(see  Figure  7).  In  these  situations,  the  downward  movement  of  water  actually  limits 
contaminant  transport  up  into  the  vadose  zone  fi'om  the  saturated  zone.  Consequently, 
since  less  mass  is  present  in  the  vadose  zone,  less  significant  recontamination  of  the 
saturated  zone  occurs.  Only  under  very  low  infiltration  rates  (<  0.04  cm/day)  does  the 
water  movement  become  too  insignificant  to  combat  the  upward  diffusion-dispersion 
contaminant  transport.  Sand  is  the  exception  which  shows  a  continual  increase  in  mass 
returned  with  infiltration  rate.  The  high  porosity  of  this  soil  allowed  diffusive-dispersive 
processes  to  dominate  throughout  the  range  of  infiltration  rates  used  in  this  case.  Under 
other  scenarios  (i.e.  ones  using  deeper  soils),  however,  the  trend  where  mass  returned 
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Figure  7:  Effect  of  Infiltration  Rate  on  Mass  Returned  to  the  Saturated  Zone 

(Vadose  Zone  Depth:  3  m) 

decreases  with  increasing  water  infiltration  is  evident.  All  these  model  restxlts  are 
included  in  Appendix  E. 

Effect  of  Vadose  Zone  Depth 

Figure  8  illustrates  the  effect  that  vadose  zone  depth  has  on  the  amount  of  mass 
returned  to  the  saturated  zone.  The  plot  clearly  shows  that  deeper  vadose  zones  cause 
more  mass  to  be  delivered  into  the  groimdwater.  This  effect  is  simply  due  to  the  fact  that 
thicker  soils  have  a  much  greater  potential  to  accumulate  contaminant  mass.  The  initial 
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Figure  8:  Effect  of  Vadose  Zone  Depth  on  Mass  Returned  to  the  Saturated  Zone 

(Infiltration  Rate:  0.20  cm/day) 


concentration  profiles  displayed  in  Appendix  C  support  this  observation  by  showing  a 
definite  increase  in  total  contaminant  mass  contained  in  the  vadose  zone  as  thicknesses 
increase.  All  the  vadose  zone  depth  results  are  contained  in  Appendix  F. 


Mass  Fraction  Returned  to  Saturated  Zone 

The  information  provided  in  the  previous  three  sections  identified  how  soil  type, 
infiltration  rate,  and  vadose  zone  depth  effected  the  amount  of  mass  returned  to  the 
saturated  zone.  However,  in  order  to  determine  what  constitutes  significant 
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recontamination,  the  mass  returned  must  be  taken  relative  to  the  mass  originally  present 

in  the  saturated  zone  (see  Equation  27): 

mass  returned  to  saturated  zone  „  .  , 

- =  mass  fraction  returned  (27) 

initial  mass  in  saturated  zone 

The  initial  masses  used  in  this  study  were  calculated  by  assuming  a  sand  aquifer  and 
using  a  reasonable  range  of  saturated  zone  thicknesses  (see  Equation  (28)). 

initial  mass  in  saturated  zone  =  Gj-Cj  (saturated  zone  thickness)  (28) 

The  mass  returned  data  was  presented  previously  in  Tables  2  and  3;  initial  masses  are 
listed  in  Table  4. 


Table  4:  Initial  Mass  Present  in  Saturated  Zone 


2.5 

Saturated  Zone  Thickness  (m) 

5.0  7.5  10.0  12.5  15.0 

Mass  (g/m^)  52\ 

10.43  16.64  20.85  26.06  31.28 

Results  illustrating  the  mass  fraction  returned  to  the  saturated  zone  are  provided  in 
Appendix  G.  For  this  study,  the  mass  fraction  threshold  was  set  at  20%.  By  using  this 
threshold,  each  scenario  was  evaluated  on  its  potential  to  significantly  recontaminate  the 
saturated  zone.  As  an  example.  Figure  9  displays  a  case  where  significant 
recontamination  is  only  likely  occur  when  a  heterogeneous  soil  {Lens  or  Layered) 
overlies  a  thin  to  medium  saturated  zone.  Figure  9  also  illustrates  that,  in  this  particular 
situation,  the  homogeneous  soils  are  unlikely  to  pose  any  problem  regardless  of  the 
thickness  of  the  saturated  zone.  All  final  results  are  categorized  in  Table  5. 


30 


Figure  9:  Mass  Fraction  Returned  to  Saturated  Zone 
(Infiltration  Rate:  0.04  cm/day;  Vadose  Zone  Depth:  3  m) 


Fortunately,  this  large  quantity  of  raw  results  can  be  reduced  to  a  few  general 
findings.  In  reference  to  causing  significant  recontamination  of  the  saturated  zone: 

1)  heterogeneous  conditions  are  worse  than  homogeneous; 

2)  the  presence  of  a  low  permeability  layer  within  a  highly  permeable  soil  seems 
to  be  the  worst  case  situation; 

3)  low  infiltration  rates  are  worse  than  high  (as  long  as  zero  is  not  approached); 

4)  deep  vadose  zones  are  worse  than  shallow;  and, 

5)  thin  saturated  zones  are  worse  than  thin 

For  a  more  quantitative  perspective,  the  reader  is  directed  to  Table  5  or  to  Appendix  G  to 
obtain  specific  information  on  a  particular  issue. 
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Table  5:  Final  Results  for  All  Scenarios 


Infil  Rate 
(cm/day) 

Sand 

Saturated  Zone  Thickness* 

Thin  Medium  Thick 

Loam 

Saturated  Zone  Thickness* 
Thin  Medium  Thick 

Vadose  Zone 
Depth 

0.00 

N 

N 

N 

N 

N 

N 

Shallow 

0.04 

N 

N 

N 

N 

N 

N 

(3  m) 

0.20 

N 

N 

N 

N 

N 

N 

0.40 

M 

N 

N 

N 

N 

N 

0.00 

S 

S 

N 

N 

N 

N 

Moderate 

0.04 

s 

s 

M 

S 

M 

N 

(10  m) 

0.20 

s 

s 

S 

s 

N 

N 

0.40 

s 

s 

s 

s 

N 

N 

0.00 

s 

s 

s 

s 

M 

N 

Deep 

0.04 

s 

s 

s 

s 

M 

N 

(30  m) 

0.20 

0.40 

s 

s 

s 

Clay* 


Infil  Rate 
(cm/day) 

Saturated  Zone  Thickness 

Thin  Medium  Thick 

Vadose  Zone 
Depth 

0.00 

S 

M 

N 

Shallow 

0.04 

S 

M 

N 

(3  m) 

0.20 

S 

M 

N 

0.40 

M 

N 

N 

Notation 

S  Significant  recontamination  (mass  factor  returned  >  threshold) 

M  Marginal  recontamination  (mass  factor  returned  o  threshold  {interval  crosses  threshold}) 
N  No  recontamination  (mass  factor  returned  <  threshold) 

-  Result  exhibited  mass  significant  balance  errors 

‘  General  thicknesses:  thin:  2.5 -5.0  m;  medium:  5.0- 10.0  m;  thick:  10.0+ m 
*  Results  for  clay  with  medium  and  deep  vadose  zones  exhibited  significant  mass  balance  errors 
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Table  5:  Final  Results  for  All  Scenarios  (cont’d) 


Lens 

Layer 

Infil  Rate 

Saturated  Zone  Thickness 

Saturated  Zone  Thickness 

Vadose  Zone 

(cm/day) 

Thin 

Medium 

Thick 

Thin 

Medimn 

Thick 

Depth 

0.00 

S 

M 

N 

S 

N 

N 

Shallow 

0.04 

S 

M 

N 

S 

N 

N 

(3  m) 

0.20 

S 

M 

N 

M 

N 

N 

0.40 

M 

N 

N 

0.00 

S 

S 

S 

S 

S 

N 

Moderate 

0.04 

S 

S 

S 

S 

M 

N 

(10  m) 

0.20 

S 

M 

N 

0.40 

S 

M 

N 

0.00 

S 

S 

s 

S 

M 

N 

Deep 

0.04 

S 

S 

s 

s 

M 

N 

(30  m) 

0.20 

0.40 

Note:  Notation  is  defined  in  the  first  part  of  this  table 


Effect  of  Model  Boundary  Conditions 

It  was  suspected  that  the  use  of  unrealistically  abrupt  lower  boundary  conditions 
in  the  model  might  have  induced  substantially  large  contaminant  fluxes  across  the  water 
table.  Since  in  the  natural  environment,  contaminant  concentrations  rise  and  fall 
gradually,  the  concern  was  that  the  model  might  significantly  overestimate  mass 
transport.  To  determine  the  extent  to  which  the  abrupt  boimdary  conditions  effected 
model  output,  four  selected  scenarios  were  ran  using  a  different  set  of  boundary 
conditions.  The  new  boundary  conditions  caused  the  aqueous  concentration  to  linearly 
rise  from  zero  to  the  specified  level  gradually  over  200  days  (vadose  zone  contamination 
phase).  Subsequently,  for  the  recontamination  phase,  the  concentration  fell  gradually  to 
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zero  over  another  200  days.  Model  results  using  these  new  gradient  boundary  conditions 
are  located  in  Appendix  H. 

Table  6  provides  a  comparison  of  the  mass  fractions  returned  under  each 
boundary  condition  scheme. 


Table  6:  Effect  of  Boimdary  Conditions  on  Mass  Fraction  Returned 
(Infiltration  Rate:  0.04  cm/day;  Vadose  Zone  Depth:  30  m) 


Mass  Fraction  Returned 

Soil  Type _ Abrupt  BC  Gradient  BC  %  Change 


Sand 

14.34 

13.47 

6.1 

Loam 

3.53 

3.43 

2.8 

Lens 

18.66 

17.83 

4.4 

Layered 

3.99 

a 

N/A 

^  "  ■■  I.-  I 

Result  exhibited  significant  mass  balance  errors 


The  small  percent  changes  shown  in  Table  6  imply  that  the  choice  of  boundary 
conditions  do  not  significantly  effect  the  model’s  outcome.  The  reason  for  this  is 
apparent  from  Figure  10  which  compares  the  fluxes  of  the  two  boundary  condition 
schemes.  Despite  greatly  different  fluxes  at  short  times,  the  overall  area  under  the  two 
flux  curves  appears  to  be  very  comparable.  As  a  result,  total  masses  determined  by 
integration  of  these  curves  would  be  similar  and  not  greatly  change  model  results. 
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Flux  Across  Water  Table  (g/m^2/day) 


Figlixe  10:  Effect  of  Boundary  Conditions  on  Flux  Across  Water  Table 

(Soil  Type:  Sand;  Infiltration  Rate;  0.04  cm/day;  Vadose  Zone  Depth:  30  m) 


This  study  had  two  purposes:  (1)  to  determine  the  fate  of  VOCs  which  migrate 
into  the  vadose  zone  from  the  saturated  zone,  and  (2)  to  determine  conditions  where 
VOCs  in  the  vadose  zone  could  cause  significant  recontamination  of  groundwater.  The 
task  was  accomplished  by  modeling  the  transport  of  TCE  between  the  saturated  and 
vadose  zones  using  differing  soil  types,  water  infiltration  rates,  and  vadose  zone  depths. 

In  general,  model  results  indicated  that  only  a  small  fraction  of  the  total 
contaminant  mass  escaped  into  the  atmosphere  or  remained  in  the  vadose  zone.  The 
majority  of  the  contaminant  mass  was  transported  back  into  the  saturated  zone.  Soil  type, 
infiltration  rate,  vadose  zone  depth,  and  saturated  zone  thickness,  all  had  a  substantial 
effect  on  the  magnitude  of  recontamination.  Deep,  heterogeneous  soil  profiles 
accompanied  with  low  infiltration  rates  posed  the  greatest  threat  to  saturated  zone 
contamination.  Moreover,  the  presence  of  a  low  permeability  lens  in  a  highly  porous  soil 
seemed  to  create  the  worst  case.  Additionally,  most  water  infiltration  rates  actually 
decreased  recontamination  of  the  groundwater  by  minimizing  the  migration  of  VOCs  into 


the  vadose  zone. 


Recommendations  for  Future  Research 

Since,  this  study  was  limited  one  dimension,  it  would  be  interesting  to  develop  a 
multi-dimensional  model  and/or  incorporate  more  realistic  boundary  conditions  to 
determine  if  the  results  significantly  change.  Efforts  to  validate  the  model  against 
existing  field  or  experimental  data  are  also  necessary. 
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Appendix  A:  Notation 


Symbol 

Description 

Dimension* 

Ca 

Solute  concentration  in  gaseous  (air)  phase 

Ci 

Initial  aqueous  phase  solute  concentration  in  saturated  zone 

Cs 

Solute  concentration  in  soil  phase  due  to  aqueous  phase  sorption 

C~^ 

Solute  concentration  in  aqueous  phase 

Da 

Effective  gaseous  phase  diffusion-dispersion  coefficient 

[L'r*] 

Dw 

Effective  aqueous  phase  diffusion-dispersion  coefficient 

[L'r‘] 

Dfa 

Free  air  molecular  diffusion  coefficient 

[L^r'] 

Dftv 

Free  water  molecular  diffusion  coefficient 

[L^r‘] 

h 

Hydraulic  head 

[L] 

hb 

Bubbling  pressure  (Brooks  and  Corey  soil  parameter) 

[L] 

J 

Solute  flux 

[M  L'^  T] 

*^srf 

Solute  flux  across  soil  surface 

[ML'^T] 

•^wt 

Solute  flux  across  water  table 

[M  L'^  T] 

K 

Unsaturated  hydraulic  conductivity 

[LT'l 

Ks 

Saturated  hydraulic  conductivity 

[Lr'l 

Kd 

SolidsAVater  distribution  coefficient 

[Lw'm;'] 

K„ 

Henry’s  law  constant  (air/water  partitioning  coefficient) 

[unitless] 

L 

Total  depth  of  vadose  zone 

[L] 

m 

Van  Genuchten  soil  parameter 

[unitless] 

n 

Effective  porosity 

[unitless] 

^wz 

Vertical  aqueous  phase  advective  flux  through  vadose  zone 

[LT’] 

*lwx 

Horizontal  aqueous  phase  advective  flux  through  saturated  zone 

[LT*] 

z 

Depth 

[L] 

a 

Van  Genuchten  soil  parameter 

[L-'] 

ttL 

Longitudinal  dispersivity  in  vadose  zone 

[L] 

tty 

Transverse  dispersivity  in  vadose  zone  (fi-om  groundwater  flow) 

[L] 

P 

Van  Genuchten  soil  parameter 

[unitless] 

0a 

Volumetric  air  content 

[unitless] 

Qr 

Residual  volumetric  water  content 

[imitless] 

0s 

Saturated  volumetric  water  content 

[unitless] 

0w 

Volumetric  water  content 

[unitless] 

X 

Pore-size  distribution  index  (Brooks  and  Corey  soil  parameter) 

[unitless] 

Pb 

Solids  bulk  density 

[MsV^] 

■ta 

Gaseous  phase  tortuosity  factor 

[unitless] 

Aqueous  phase  tortuosity  factor 

[unitless] 

'F 

Pressure  head 

[L] 

♦  Dimensions;  L  =  length 

Dimension  subscripts:  a  =  air 

V  =  voids 

M  =  mass 

m  =  media 

w  =  water 

T  =  time 

s  =  solids 
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0.50 


Water  Content  vs.  Depth 
Loam  (Vadose  Zone  Depth  =  3  m) 
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■2 


Water  Content  vs.  Depth 

Layered  Soil  (Vadose  Zone  Depth  =  3  m) 
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Water  Content  vs.  Depth 
Layered  Soil  (Vadose  Zone  Depth  =  30  m) 


Infiltration 
Rate  (cm/day) 


10  15  20 

Depth  (m) 


Water  Content  vs.  Depth 

Layered  Soil  (Vadose  Zone  Depth  =  30  m) 
0.50  - 

0.40  - 


0.10 


0.00  T-h . » . . 4  >■'■■■  > . I . t . t . >■■■-< . +-+-f . I . I  -■I . t . t . < . 


Infiltration 
Rate  (cm/day) 


- 0.04 


0  5  10  15  20  25  30 

Depth  (m) 


0.50 
0.40 ; 

s 

|o.3o ; 

o 

o 

|o.20  ; 
0.10  ; 


Water  Content  vs.  Depth 
Layered  Soil  (Vadose  Zone  Depth  =  30  m) 


Infiltration 
Rate  (cm/day) 


0.00 


Water  Content  vs.  Depth 

Layered  Soil  (Vadose  Zone  Depth  =  30  m) 


Infiltration 
Rate  (cm/day) 


Appendix  C:  Initial  Vadose  Zone  Concentration  Profiles 


Relative  Total  Concentration  vs.  Depth 


Sand  (Vadose  Zone  Depth  =  3  m) 


Depth  (m) 


Infiltration 
Rate  (cm/day) 


- 0.00 

- -.0.04 

—  ■  -  —  0.20 
- 0.40 


Relative  Total  Concentration  vs.  Depth 
Sand  (Vadose  Zone  Depth  =  10  m) 


Infiltration 
Rate  (cm/day) 


- 0.00 

. 0.04 

-  ■ "  -  0.20 
— - 0.40 


Depth  (m) 


Relative  Total  Concentration  vs.  Depth 
Sand  (Vadose  Zone  Depth  =  30  m) 
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Relative  Total  Concentration  vs.  Depth 

Loam  (Vadose  Zone  Depth  =  3  m) 


l.OE+00 
l.OE-01 

g  l.OE-02 
l.OE-03 

3 

j2  l.OE-04 

>  l.OE-05 

«3 

^  l.OE-06 
l.OE-07 
l.OE-08 

0.0  0.5  1.0  1.5  2.0  2.5  3.0 

Depth  (m) 


Relative  Total  Concentration  vs.  Depth 

Loam  (Vadose  Zone  Depth  =  10  m) 

l.OE+OO 
l.OE-01 
^  l.OE-02 
J  l.OE-03 
3  l.OE-04 
^  l.OE-05 
I  l.OE-06 
&  l.OE-07 
l.OE-08 
l.OE-09 

0  2  4  6  8  10 

Depth  (m) 


Relative  Total  Concentration  vs.  Depth 

Loam  (Vadose  Zone  Depth  =  30  m) 


e 

5 

3 

f2 

u 

.> 


l.OE+OO 
l.OE-02 
l.OE-04  T 
LOE-06 
l.OE-08  + 
l.OE-10 
l.OE-12  + 
l.OE-14 
l.OE-16 
l.OE-18  ' 

LOE-20 


Infiltration 
Rate  (cm/day) 

-0.00 
-0.04 

0.20  and  0.40 
are  not  shown 
due  to  errors  in 
mass  balance 


15  20 

Depth  (m) 


Relative  Total  Concentration  vs.  Depth 
Clay  (Vadose  Zone  Depth  =  3  m) 


Not  shown  are  the  concentration  profiles  for  clay  with  vadose 
zone  depths  of  10  and  30  m.  Those  cases  exhibited  significant 
mass  balance  errors. 
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Relative  Total  Concentration  vs.  Depth 
Sand  w/  0.5  m  Clay  Lens  (Vadose  Zone  Depth  =  3  m) 


Depth  (m) 


Infiltration 
Rate  (cra/day) 

■■-0.00 

. 0.04 

- 0.20 


♦  0.40  is  not 
shown  due  to 
errors  in  mass 
balance 


Relative  Total  Concentration  vs.  Depth 
Sand  w/  0.5  m  Clay  Lens  (Vadose  Zone  Depth  =  10  m) 


Relative  Total  Concentration  vs.  Depth 

Sand  w/  0.5  m  Clay  Lens  (Vadose  Zone  Depth  =  30  m) 
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Relative  Total  Cone.  Relative  Total  Cone.  Relative  Total  Cone 


Mass  (g/m''2)  Mass  (g/m''2) 


Fate  of  Contaminant  Mass  vs.  Soil  Type 
(Infiltration  Rate  =  0.20  cm/day;  Vadose  Zone  Depth  =  3  m) 


Not  shown  are  the  results  with  an  infiltration  rate  of 
0.20  cm/day  and  a  vadose  zone  depth  of  30.  For  this  case, 
all  soils  exhibited  significant  mass  balance  errors. 


Fate  of  Contaminant  Mass  vs.  Soil  Type 
(Infiltration  Rate  =  0.40  cm/day;  Vadose  Zone  Depth  =  3  m) 


i 


□Atmosphere 
□Vadose  Zone 
■  Sat.  Zone 


♦  Data  for  clay 
and  lens  are  not 
shown  due  to 
errors  in  mass 
balance 


Sand  Loam  Clay*  Lens*  Layered 
Soil  Type 


Fate  of  Contaminant  Mass  vs.  Soil  Type 

(Infiltration  Rate  =  0.40  cm/day;  Vadose  Zone  Depth  =  10  m) 


□Atmosphere 
□Vadose  Zone 
■  Sat.  Zone 


*  Data  for  clay 
and  lens  are  not 
shown  due  to 
errors  in  mass 
balance 

Sand  Loam  Clay*  Lens*  Layered 
Soil  Type 


1 


8.0 

7.0 

6.0 

5.0 

4.0 

3.0 

2.0 

1.0 

0.0 


Fate  of  Contaminant  Mass  vs.  Soil  Type 

(Infiltration  Rate  =  0.40  cm/day;  Vadose  Zone  Depth  =  30  m) 


□Atmosphere 
□Vadose  Zone 
■  Sat.  Zone 


*  Data  for  these 
soils  are  not 
shown  due  to 
errors  in  mass 
balance 


Mass  Returned  to  Saturated  Zone  vs.  Infiltration  Rate 


— • — 

-Sand 

-Loam 

- IK— 

-Clay 

--tr- 

■Lens 

— e— 

-Layered 

0.0  0.1  0.2  0.3  0.4 

Infiltration  Rate  (cm/day) 


Mass  Returned  to  Saturated  Zone  vs.  Infiltration  Rate 


0.0  0.1  0.2  0.3  0.4 

Infiltration  Rate  (cm/day) 


-Sand 

■  -  *  - 

•Loam 

■  ■  A-  • 

•Lens* 

— ©— 

-Layered 

*  Data  for  clay 
and  lens  are  not 
shown  due  to 
errors  in  mass 
balance 


Mass  Returned  to  Saturated  Zone  vs.  Infiltration  Rate 


(Vadose  Zone  Depth  =  30  m) 


—4 — 

-Sand* 

•Loam* 

-  “  -A-  “ 

•Lens* 

— ©_ 

-Layered* 

*  Some  data  are 
not  shown  due 
to  errors  in 
mass  balance 
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Mass  Returned  to  Sat.  Zone  vs.  Vadose  Zone  Depth 
(Infiltration  Rate  =  0.00  cm/day) 


Appendix  G:  Mass  Fraction  Returned  to  the  Saturated  Zone 


Fraction  of  Mass  Returned  to  Mass  Initially  Present  in  Sat 
Zone  vs.  Sat  Zone  Depth 
(Infiltration  Rate  =  0.00  em/day/Vadose  Zone  Depth  =  3  m) 


I 

ed 


n 

s 


Fraction  of  Mass  Returned  to  Mass  Initially  Present  in  Sat. 
Zone  vs.  Sat  Zone  Depth 
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Fraction  of  Mass  Returned  to  Mass  Initially  Present  in  Sat. 
Zone  vs.  Sat.  Zone  Depth 
(Infiltration  Rate  =  0.04  cm/dayA^adose  Zone  Depth  =  3  m) 
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Fraction  of  Mass  Returned  to  Mass  Initially  Present  in  Sat. 
Zone  vs.  Sat.  Zone  Depth 
(Infiltration  Rate  =  0.20  cm/day/V adose  Zone  Depth  =  3  m) 
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Saturated  Zone  Thickness  (m) 


Not  shown  are  the  results  with  an  infiltration  rate  of 
0.20  cm/day  and  a  vadose  zone  depth  of  30.  For  this  case, 
all  soils  exhibited  significant  mass  balance  errors. 


Fraction  of  Mass  Returned  to  Mass  Initially  Present  in  Sat 
Zone  vs.  Sat  Zone  Depth 
(Infiltration  Rate  =  0.40  cm/dayA'adose  Zone  Depth  =  3  m) 
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Fraction  of  Mass  Returned  to  Mass  Initially  Present  in  Sat. 
Zone  vs.  Sat.  Zone  Depth 
(Infiltration  Rate  =  0.40  cm/dayA^adose  Zone  Depth  =  10  m) 


Fraction  of  Mass  Returned  to  Mass  Initially  Present  in  Sat. 
Zone  vs.  Sat.  Zone  Depth 
(Infiltration  Rate  =  0.40  cm/dayA^adose  Zone  Depth  =  30  m) 


♦  Data  for  all 
other  soils  are 
not  shown  due 
to  errors  in 
mass  balance 
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Initial  Vadose  Zone  Concentration  Profiles: 


Relative  Total  Concentration  vs.  Depth 

(Infiltration  Rate  =  0.04  cm/day;  Vadose  Zone  Depth  =  30  m) 

l.OE+00 
g  l.OE-01 

O 

U 

3 

j2  l.OE-02 

u 
to 

3  l.OE-03 
cc 

l.OE-04 

0  5  10  15  20  25  30 

E>epth  (m) 


Relative  Total  Concentration  vs.  Depth 
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not  shown  due 
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mass  balance 


Not  shown  are  the  concentration  profiles  for  the  clay  and 
layered  cases  which  exhibited  significant  mass  balance  errors. 
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Mass  (g/m'^2) 


Appendix  I:  Equation  Development  for  Soil  Water  Content  Model 


1.  Apply  Buckingham  Equation  to  define  'E(z): 

a.  Assumptions: 

1.  Discontinuous  porous  medium  can  be  described  by  a  hypothetical  continuous 
representative  elementary  volume  (REV) 

2.  One  dimensional  laminar  flow 

3.  Porous  media  is  homogeneous,  isotropic,  isothermic,  and  nondeformable 

4.  Fluids  are  incompressible  (densities  of  water  and  air  are  constant) 

5.  Steady  state 

6.  Velocity  head  is  negligible 

b.  Equation  Development: 

q^  =  -K('P)—  where:  h  =  ^  +  z  (Buckingham,  1907) 

3z 


c.  Finite-difference  approximation  to  be  solved  interatively  for  'P(z): 


where: 


'I'avg  2 

K  Jl  -  (a'P„,)'>-‘[l  + 

K('Pavg)  =  — - - m - -  Genuchten,  1980) 

““(KT'  P  =  >.+i 

d.  Boundary  Condition:  'P(z  =  L)  =  0 
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2.  Solve  for  0w(z): 


0^(z)  =  0r  + 


i  +  [aT(z)]P]'" 


(van  Genuchten,  1980) 
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Appendix  J:  Equation  Development  for  Contaminant  Transport  Model 


L  Fundamental  Mass  Transport  and  Flux  Equations  for  the  Unsaturated  Zone: 

a.  Assumptions: 

1 .  Discontinuous  porous  medium  can  be  described  by  a  hypothetical  continuous 
representative  elementary  volume  (REV) 

2.  One  dimensional  laminar  flow 

3.  Fluids  are  incompressible  (densities  of  water  and  air  are  constant) 

4.  No  NAPL  phase  present 

5.  Vapor  sorption  is  negligible 

6.  Air  phase  is  passive  (gaseous  advective  flux  is  negligible) 

7.  All  other  sources  and  sinks  are  negligible 

b.  Contaminant  Transport  Equation; 

c.  Flux  Equation: 

2.  Simplify  Equation: 
a.  Assumptions: 

1 .  Porous  media  is  homogeneous,  isotropic,  and  isothermic  (Df^,  Df^^  and  n  are 
constant  in  space) 

2.  Porous  media  is  nondeformable  (pj,  and  n  are  constant  in  time) 

3.  Steady  state  vertical  water  flow  (q^^  and  is  constant  with  time  and  space;  0^  and 
0a  are  constant  in  time) 

4.  Local  equilibrium 

5.  Solids/water  distribution  and  air/water  partitioning  are  linear: 

Cs  = 

Ca  =  KjjC^  (assuming  ideal  behavior  =>  aqueous  activity  coefficient  =1) 

6.  Millington  formula  is  appropriate  to  estimate  and  Xa; 

7  7 

Xw  =  ——2^  Xa  =  ^^2  (Millington,  1959) 

n  n 
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b. 


Contaminant  Transport  Equation: 


ac^  _ 

f'  \ 

qwz 

,  1  5 

at 

IRwJ 

dz  Ry,  dz 

(D„+D.K„)^ 


where: 


Rw  =0w+0aKH+PbKd 


Dw  = 


Dfw 


Mwx 


+  a, 


ew^ 


le^T^  and:  D,  =  DfA'^a 


0a  =  n-0^ 


=  K(0^)—  ( —  is  the  hydraulic  gradient 

dx.  V5x 


c.  Flux  Equation: 

J  =  qwz^w  ~  (Dw  +  DjKh) 


dz 


3.  Finite-Difference  Numerical  Approximation  of  Contaminant  Transport  Equation 
(Crank-Nicolson  Method  for  Diffusion/Central  Space  Method  for  Advection): 


~ + — + D,K„)(c,;'.,  -  C„;' + -  Cm' 

2Rwi(Az^ 

D.K„),_,(c„i'  +c„"' 


1  n 

j 

n  -1- 

^wi-1  "^^Wi+I  J 

'wi 

4R^jAz 

M _ 

At 

T - 

2R 

wi(Azr 

-(D„  +D.K„),.,(c„r  +C."‘  -C„"’)] 

at 


R 


Wl- 


-(Dw- 
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d.  Simplify  Equation: 

T  -V  At 

Let:  A,,  = - 

h  4R^jAz 


and;  A,.  = 


At 


2R„,(Azy 


and:  D,=(D„+D.K„) 


Urns:  C„”*' 


—  o  ®— n  1  tf'  **  _  p  **  4.  r*  _  r’ 

“  '^Wi  4wz^l.  Vrwi+1  Wi-1  '^wi+i  Wi-1  J 

-D,.A,.(c„f-c„i'_,+c„r‘-c„":,') 


e.  Move  All  “n  +1”  Terms  to  Left  Side  of  Equation; 

-(q^X,.  +[l  +  X^^(D,  +D,_,)]c„«' 

=  +  [i  -  (D, + D,_,)]c.j  -  (q„x,,  -  D,x^ 


n+1 

i+1 


f.  Rewrite  Solution  as  a  Tridiagonal  System  of  Equations: 


bi 

ei 

0 

0 

. 

0 

0 

Cwj 

RHSi  -  ai  Cwo"^* 

a2 

b2 
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0 

0 

• 

0 

Cw2 

RHS2 

0 

as 

^3 
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. 

» 

Cw3 

RHS3 

0 

0 

• 

” 

. 

• 

« 

. 

. 

. 

.  .  . 

0 

0 

. 

• 

• 

• 

• 

^ta-3  ^nz-3 

6nz-3 

0 

Cwnz.3 

Rhs^.3 

0 

» 

• 

« 

^nz-2 

6nz-2 

Cw„z.2 

RHSnz.2 

0 

0 

• 

• 

•  0  0 

anz-1 

bnz-l 

Cw„z.l 

RHS„,.i-e^.iCwJ^' 

where: 


1  +  A,  (Dj  +Dj_i) 


RHS,=-a,C„;.,  +  ^-b,):„l' 


-e|C 


n 

wi+i 


Note:  The  system  is  solved  using  the  Thomas  algorithm.  The  condition, 
t>i  >  N I  +  Icj  I  >  0,  was  used  to  ensure  that  the  coefficient  matrix  was 
diagonally  dominant  (Mayers  and  Morton,  1994). 
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4.  Finite-Difference  Numerical  Approximations  of  Flux  Equation  for  Fluxes  at  Surface 
and  at  Water  Table  (Backward  Space  Method): 

Jsrf=q»zC„(z  =  0)  -  5^^  (C,  (z  =  0)  -  C„  (z  =  0  -  Az)) 

Az 

=  q™C„  (z  =  L)  -  5^^(C„(z  =  L  -  Az)  -  C„  (z  =  L) 

Az 


5.  Inital  and  Boundary  Conditions: 

IC:  Cw(z,  t  =  0)  =  0  BCl:  Cw(z  =  0,t)  =  0  {at  surface} 

BC2:  C^(z  =  L,  t  =  0..tf,*)  =  Ci  {at  water  table} 
Cw(z  =  L,  t  =  tfi ..tf2*)  =  0  {at  water  table} 

*tfi  is  2000  days 

tfj  is  the  time  it  takes  for  the  highest  aqueous  concentration  in  the  system  to  fall 
below  1  pg/L  or  2000  days 
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Appendix  K:  Analytical  Solution  to  Simplified  Contaminant  Transport  Equation 


The  following  derivation  was  obtained  from  a  personal  communication  with  D.  Quinn, 
Department  of  Mathematics  and  Statistics,  Air  Force  Institute  of  Technology,  OH, 

18  November  1996.  The  method  mirrors  a  dififiision-only  approach  described  in 
Haberman  (1987). 


1 .  Simplify  mass  transport  equation  by  assuming  coefficients  for  advection,  U,  and 
diffusion-dispersion,  D,  are  constant: 


general  equation:  =  -U  +  D 

with  inital  condition:  C^(z,0)  =  0 
and  boimdary  conditions:  C^(0,t)  =  0 


dz^ 


Cw(L,t)  =  C, 


2.  Transform  C^(z,t)  to  produce  homogeneous  boimdary  conditions: 


define:  (z,  t)  =  —  z  +  V (z,  t) 

L 


5V  C-  5V 

so  equation  becomes:  —  =  -U  —  -  U  —  +  D  — ^ 

at  L  az  az^ 


C- 

with  inital  condition:  V(z,0)  =  — ^^z 


and  boundary  conditions:  V(0,t)  =  0  V(L,t)  =  0 


av 

3.  Transform  V(z,t)  to  eliminate  the  —  term  using  a  form  suitable  for  a  Fourier  series: 

az 


define:  V(z,t)  =  e2DW(z,t) 
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.  ^  aw 

so  equation  becomes:  - 

at 


UC; 

- *-e 

L 


ife 

2D  _ 


4D 


W  +  D 


a^w 

az^ 


with  inital  condition:  W(z,0)  =  - 


Uz 

Cz 

2d 


and  boundary  conditions:  W(0,t)  =  0  W(L,t)  =  0 


4.  Represent  W(z,t)  as  a  Fourier  series  which  satisfies  the  boundary  conditions 
W(0,t)  =  0andW(L,t)  =  0: 

W(z,t)  =  2]a„(t)sin-— 

n=l  ^ 


SO  equation  becomes: 

nTiz 


S“n'(t) 
n=l 


sin- 


UCi  .  nTTz 

L  4Dft  "  L 


n=l  ^ 


2  2 

n  7t  .  niiz 
:5— sin - 

2  T 


kitz 

5.  Multiply  through  transport  equation  by  sin - and  integrate  from  0  to  L: 

L 


Z.  .  ^  f  ,  nTtz  ,  kirz , 

“n  (t)Jsin-— sin— dz  = 


UCj  .  kiiz  A  .  niiz  .  kirz , 

- e  -^*^sin— dz - >  a_  sm - sm— dz 

L  J  L  "J  L  L 

n^Tt^  V  .  miz  .  kitz . 

-  Jsin— sm-^z 

n=l  ^  0  ^  ^ 


,  f  .  mtz  .  kTTZ .  —  if  n  =  k 

where:  sm - sm— dz=2 

J  T  T 

0  ^  ^  0  if  n^tk 


so  equation  becomes:  aj^'  (t)  =  - 


2UCi 


L  Uz 

Je"2D 

0 


.  kiiz , 
sm— dz-tti, 
L  ^ 


kVl 
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8.  Solution; 

C^(z,t)  =  ^z  +  V(z,t) 

JU 

where: 


V(z,t)  =  e2DW(z,t) 


W(z,t)  =  2ak(t)sin^ 


k=l 


a^(t)=a^(0)e-“'‘‘+It(l-e-“‘') 

Wt  ''  ^ 


...  2k7iCi 

«k(0)  =  -rr^ 


u 

2D  DL 


(  UL 


UL 

e  cos(k7i) 
k^Tl^ 


e  2®cos(k7t)-l 


V 


4D^  \} 


kV' 

+  ■ 


SkTtP^UCj 
~  4kVD2L  +  U^L^ 


/'  UL 


e  ^^cos(k7r)-l 


2„2 


cOk  =  —  +  D — ^ 
^  4D 
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Appendix  L:  FORTRAN  computer  code  for  Soil  Water  Content  Model 


************************************************************************ 

*  PROGRAM:  MODELWC  -  Determines  Pressure  Head,  Conductivity,  and  * 

*  Water  Content  as  a  function  of  Depth  * 

*  -  Output  file  contains  POR{I)  and  PB(I)  which  is  * 

*  needed  in  MODELCT  {note:  WCS=P0R}  * 

*  -  PSI  is  defined  as  negative  pressure  head,  so  * 

*  PSI (I)  is  positive  * 

************************************************************************ 

C 

C 

C  Variable  Declarations 

************************************************************************ 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z) 

DIMENSION  PSI (0:1000) ,  ZK{0:1000),  WC(0:1000) 

CHARACTER*10  INFIL,  OUTFIL,  INPUT,  OUTPUT 

DOUBLE  PRECISION  PSIA,TZK,TQZ,DQZO,DQZN,TOL,DELZ,VA,VB,VM, 

@QZ , WCS, WCR, ZKS, HB, CLAM, PB, ZL, T, Z 
COMMON  QZ , WCS , WCR, ZKS , HB , CLAM, PB, ZL, T, Z 
C 
C 

C  Format  Statements 

************************************************************************ 

C 

1000  FORMAT (AlO) 

2000  FORMAT (' MODEL  WC:  Pressure  head  and  water  content  vs  depth',//, 

® 'OUTPUT: ',/, 3X, ' Z=Depth;  PSI=Pressure  Head;  K=Conductivity; 
@WC=Water  Content;  POR=Porosity;  PB=Soil  Bulk  Density',//, 

@7X, 'Z' ,10X, 'PSI' ,10X, 'K' ,10X, 'WC , 9X, 'POR' ,10X, 'PB') 

3000  FORMAT (6E12 , 5) 

4000  FORMAT (IX//, ' Solution  Does  Not  Converge--Try  Larger  nz ' ) 

C 


Define  Input/Output  Files 

icieie^’kicitiritit^it’k'kieie’k'k'k’k’k'k'k'k'k'k'kic'kie'ie'k'kie'kie’k'k'k'k'k'kic'kieifif'kieie'k'k'k'k-kieie'kie'k'k'kieie'k'k'k'k'kiciif’k 

c 

WRITE (6,*)  'Name  of  input  file' 

READ (5, 1000)  INFIL 
INPUT  =  INFIL 

WRITE (6,*)  'Name  of  output  file' 

READ (5, 1000)  OUTFIL 
OUTPUT  =  OUTFIL 

OPEN  (UNIT  =  10,  FILE  =  OUTPUT,  STATUS  =  'NEW') 

OPEN  (UNIT  =  20,  FILE  =  INPUT,  STATUS  =  'OLD') 

WRITE  (10,  2000) 
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on  n  o  o  no 


C 

C 

C  Calculations 

************************************************************************ 

C 

C  Initialize  variables 
C 

TOL=0. 000001 

MAXIT=1000 

NZ=0 

PSI (NZ) =0.d0 
READ (20,*)  NSL,ZL,QZ 
DO  200  K=1,NSL 
C 

C  Define  Boundary  Condition;  Layer  Above  Water  Table:  PSI{Z=L)=0; 

C  Subsequent  Layers ; 

PSInew(Z=Z0)=PSIold(Z=Znz) 


PSI(0)=PSI(NZ) 

READ (20, *)  WCS,WCR, ZKS , HB, CLAM, PB, T, NZ 

Define  Parameters 

DELZ=T/NZ 
VA=l.d0/HB 
VB=CLAM+l.d0 
VM=CLAM/VB 
IF(K.EQ.l)  Z=ZL 
IF(K.EQ.l)  ZK(0)=ZKS 
IF(K.EQ.l)  WC(0)=WCS 

IF(K.EQ.l)  WRITEdO,  3000)  Z,  PSI(O),  ZK(0),  WC(0),  WCS,  PB 
IF(K.EQ.l)  PRINT*,  Z,  PSI (0) ,  WC(0) 

For  seed  value,  set  PSIA=PSI (I-l)  and  calculate  TZK 

DO  50  1=1, NZ 
COUNT=0 
PSIA=PSI (I-l) 

TZK=ZKS* ( (l.dO- (VA*PSIA)** (VB-l.dO) * (l.d0+ (VA*PSIA) 
@**VB) ** (-VM) ) **2.d0) / ( (l.d0+ (VA*PSIA) **VB) ** (VM/2 .dO) ) 
GOTO  10 
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c 

C  Iteration  Loop 

C  Note;  Since  the  down  direction  is  defined  to  be  positive,  but  the 
C  model  calculates  from  the  water  table  up  (-z  direction)  s 

C  1.  DELZ  is  negative  (with  DELZ  Z  defined  as  Z(1)-Z(I-1)) 

C  2,  TQZ  is  negative  (to  calculate  a  positive  TQZ,  the  sign 

C  is  changed  for  the  right  hand  sign  of  the  equation) 

C 

C  The  above  comments  are  reflected  in  the  equations  for  PSI(I) 

C  and  TQZ  below: 

C 

10  COUNT=COUNT+l 

IF (COUNT. GT.MAXIT)  GOTO  100 

PSI (I) = (qz/TZK-l.dO) * (-DELZ)+PSI (I-l) 

IF(PSI(I) .LT.O)  GOTO  100 
PSIA=(PSI(I)+PSI(I-1) )/2.d0 

TZK=ZKS* ( (l.dO- (VA*PSIA) ** (VB-l.dO) * (l,d0+ (VA*PSIA) 

®**VB) ** (-VM) ) **2 .do) / ( (1 .d0+ (VA*PSIA) **VB) ** (VM/2 .dO) ) 

TQZ=TZK*( (PSI (I) -PSI (I-l) )/(-DELZ)+l.d0) 

DQZN=QZ-TQZ 

IF(DQZN.LT.O)  DQZN=-DQZN 
DQZO=DQZN 

IF(DQZO.GT.TOL)  GOTO  10 
C 

C  Calculate  K  and  WC,  and  Print  Results 
C 

ZK(I) =TZK 

WC (I) =WCR+ (WCS-WCR) / ( (l.d0+ (VA* (PSI (I) ) ) **VB) **VM) 

Z=Z-DELZ 

WRITEdO,  3000)  Z,  PSI(I),  ZK(I),  WC(I),  WCS,  PB 
PRINT*,  Z,  PSI (I),  WC(I) 

50  CONTINUE 

GOTO  200 
C 

100  WRITEdO,  4000) 

PRINT*,  'Solution  Does  Not  Converge- -Try  Larger  nz ' 

GOTO  500 
C 

200  CONTINUE 
GOTO  500 
C 

500  END 
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non  non 


C 

C  Iteration  Loop 

C  Note:  Since  the  down  direction  is  defined  to  be  positive,  but  the 
C  model  calculates  from  the  water  table  up  (-z  direction) ; 

C  1.  DELZ  is  negative  (with  DELZ  Z  defined  as  Z(I)-Z(I-1)) 

C  2.  TQZ  is  negative  (to  calculate  a  positive  TQZ,  the  sign 

C  is  changed  for  the  right  hand  sign  of  the  equation) 

C 

The  above  comments  are  reflected  in  the  equations  for  PSI(I) 
and  TQZ  below: 

10  COUNT=COUNT+l 

IF (COUNT. GT.MAXIT)  GOTO  100 

PSI (I) = (qz/TZK-l.dO) * (-DELZ) +PSI (I-l) 

IF(PSI(I) .LT.O)  GOTO  100 
PSIA= (PSI (I) +PSI (I-l) ) /2.d0 

TZK=ZKS* ( (l.dO- (VA*PSIA) ** (VB-l.dO) * (1 . d0+ (VA*PS1A) 

@**VB) ** (-VM) ) **2 .do) / ( (1 .d0+ (VA*PSIA) **VB) ** (VM/2 .dO) ) 

TQZ=TZK* ( (PSI (I) -PSI (I-l) ) / (-DELZ) +l.d0) 

DQZN=QZ-TQZ 

IF(DQZN.LT.O)  DQZN=-DQZN 
DQZO=DQZN 

IP(DQZO.GT.TOL)  GOTO  10 
Calculate  K  and  WC,  and  Print  Results 
ZK(I) =TZK 

WC(I)=WCR+(WCS-WCR)/( (l.d0+(VA*(PSI(I) ) ) **VB) **YM) 

Z=Z-DELZ 

WRITE (10,  3000)  Z,  PSI (I),  ZK(I) ,  WC(I),  WCS,  PB 
PRINT*,  Z,  PSI (I),  WC(I) 

50  CONTINUE 

GOTO  200 

I 

100  WRITE (10,  4000) 

PRINT*,  'Solution  Does  Not  Converge- -Try  Larger  nz ' 

GOTO  500 

200  CONTINUE 
GOTO  500 

500  END 
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Appendix  M:  FORTRAN  computer  code  for  Contaminant  Transport  Model 


************************************^*********************************** 


*  PROGRAM;  MODELCT 

* 

* 

* 

* 

* 

* 

4 


1-D  ADVECTIVE-DISPERSIVE  TRANSPORT  IN  VADOSE  ZONE  * 
Crank-Nicolson  approximation  for  concentration  * 
Backward  space  approximation  for  flux  * 
Z,  K,  WC,  Pb,  and  porosity  are  read  in  from  an  * 
input  file  generated  from  WC  model)  * 
Since  system  is  stiff,  allows  time  domain  to  be  * 
divided  up  into  different  regions  with  different  * 
DEBT  * 


*  -  Calculates  Cw,  Ct,  J(z=0),  and  J(z=L)  * 

************************************************************************ 


*  Note  to  the  reader:  This  code  actually  contains  two  models,  one  with* 

*  fixed  boundary  conditions  (BCs)  and  one  with  * 

*  gradient  BCs.  For  the  fixed  BC  model,  activate  * 

*  the  Cl  lines,  for  the  gradient  model,  activae  * 

*  activate  the  C2  lines.  * 

****★*****★***★★★*****★*****★★***★*★****★****★*****★*★****************★* 


C 

C 

C  Variable  Declarations 

************************************************************************ 


C 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z) 

DIMENSION  CO (0:450) ,  C{0:450),  Z(0;450),  WC(0:450),  ZK{0:450), 
@POR(0:450),  PB(0:450),  RW(0:450),  D(0:450),  A(0:450),  B{0:450), 
®E(0:450),  RHS(0:450),  G(0;450),  P(0:450),  R(0;450) 

CHARACTER*10  INFILl,  INFIL2 ,  OUTFIL,  INPUTl,  INPXjr2 ,  OUTPUT 
DOUBLE  PRECISION  DELT, TO, AC, TW, TA,QWX, CLAMl , CLAM2 , SA, SE, STAB, 

®CT , JSF , JWT , Cl , GRAD , DFW , DFA, ZKD , ZKH , AT , AL , QWZ , Z 
COMMON  C I , GRAD , DFW , DFA , ZKD , ZKH , AT , AL , QWZ 

C 

C 

C  Format  Statements 

if*'k’kicicic*icieieititie*ieifieirie*iririr*ieic'k'kieit'k**'k'kir'k'k’k*'k’k’k'k-k’kic'k'k'k-k-k-k'kiric'kiei('ir'k'kici(i['kie'kiric’k 

c 

1000  FORMAT (AlO) 

2000  FORMAT (IX 'MODEL  CT:  1-D  Transport  Through  Vadose  Zone') 

3000  FORMAT (IX//, 'LOADING  OUTPUT;') 

3100  FORMAT (IX//, 'UNLOADING  OUTPUT;') 

3500  FORMAT (IX/, ' Stability  Requirements :’,/, IX, ' b- (] a | + | c ]) >0 ' 

®,1X, ' la|+|c|>0') 

3600  FORMAT (2E12. 5) 

4000  FORMAT ( IX/, 'T=Time;  Z=Depth;  Cw=Water  Cone;  Ct=Total  Cone;',/, 
®'Jsf=Flux  at  Surface;  Jwt=Flux  at  Water  Table',//, 

®7X, 'T' ,11X, 'Z' ,10X, 'Cw' ,10X, 'Ct' ,10X, ' Jsf ' ,9X, ' Jwt' ) 

5000  FORMAT (6E12. 5) 
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o  o 


C  Define  Input /Output  Files 
C 

WRITE (6,*)  'Name  of  Parameter  input  file' 

READ (5, 1000)  INFILl 
INPUTl  =  INFILl 

WRITE (6,*)  'Name  of  Water  Content  input  file' 
READ (5, 1000)  INFIL2 
INPUT2  =  INFIL2 

WRITE (6,*)  'Name  of  output  file' 

READ (5, 1000)  OUTFIL 
OUTPUT  =  OUTFIL 


OPEN 

(UNIT  =  10, 

FILE  =  OUTPUT, 

STATUS  = 

'NEW 

) 

OPEN 

(UNIT  =  20, 

FILE  =  INPUTl, 

STATUS  = 

'OLD 

) 

OPEN 

(tJNIT  =  25, 

FILE  =  INPUT2, 

STATUS  = 

'OLD 

) 

WRITE (10, 2000) 

READ (20,*)  CI,GRAD,DFW,DFA, ZKD,ZKH 
READ(20,*)  NZ,NTR,NTL,AT,AL,QWZ 
DO  10  1=0, NZ 

READ  (25,*)  Z(I) ,ZK(I) ,WC(I) ,POR(I) ,PB(I) 

10  CONTINUE 
C 

C  Define  Parameters 
C 

DO  20  1=0, NZ 

AC=POR(I) -WC(I) 

TW= (WC(I) ** (7,d0/3.d0) ) / (POR(I) **2.d0) 

TA= (AC** (7.d0/3 ,d0) ) / (POR(I) **2 .dO) 

QWX=ZK(I) *GRAD 

D (I) = (DFW+AT*QWX/WC (I) +AL*QWZ/WC (I) ) *WC (I) *TW+DFA*AC*TA*ZKH 
RW (I) =WC (I) +AC*ZKH+PB (I) *ZKD 
20  CONTINUE 
C 

C  Begin  Time  Region  Loop 
C 

T0=0 

DO  900  U=1,NTR 
IF(U.LE.NTL)  WRITE (10, 3000) 

IF(U.GT.NTL)  WRITE (10, 3100) 

WRITE (10, 3500) 

READ  (20,*)  TF,NT 
DELT=TF/NT 
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C 

C  Define  Parameters  which  are  Dependent  on  DELT 
C 

DO  50  I=1,NZ-1 

CLAM1=DELT/ (4 .dO*RW(I) * (Z (I) -Z (I-l) ) ) 
CLAM2=DELT/ (2 .dO*RW (I) * ( (Z (I) -Z (I-l) ) **2 .dO) ) 
A(I)=-  (QWZ*CIiAMl+D{I-l)  *CLAM2) 
B(I)=l.dO+CLAM2*  (Dd)  +D(I-1) ) 

E (I) =QWZ*CLAM1-D (I) *CLAM2 
IF(A(I) .LT.O)  SA=-A(I) 

IF{A(I) .GE.O)  SA=A(I) 

IF(E(I) .LT.O)  SE=-E(I) 

IF(E(I) .GE.O)  SE=E(I) 

STAB=B(I) - (SA+SE) 

WRITE (10, 3600)  STAB,  SA+SE 
50  CONTINUE 

WRITE (10, 4000) 

IF(U.EQ.l)  GOTO  100 
IF(U.EQ.NTL+1)  GOTO  200 
GOTO  500 


C  Define  Initial  Condition 

★★it********************************************************************* 

c 

C  LOADING  Phase:  C(Z,T=0)=0 

C  Note:  Since  CO  is  zero,  no  iinit  conversion  is  necessary 
C 

100  T=0.d0 

CT=0.d0 
JSF=0.d0 
JWT=0.d0 
DO  120  1=0, NZ 
CO(I) =0.d0 

WRITE (10,  5000)  T,  Z(I),  CO(I),  CT,  JSF,  JWT 
PRINT*,  T,  Z(I),  CO (I) 

120  CONTINUE 

GOTO  300 
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c 

C  UNLOADING  Phase;  C(Z,T=TF1) 

C  {final  result  from  LOADING  Phase  except  C (Z=L, T=TF1) =0} 

C  Note;  Convert  CO  back  to  ppm  and  J  to  g/m*2  for  printing 

C 

Cl  200  C0(0)=0.d0 

Cl  DO  220  1=0, NZ 

C2  200  DO  220  1=0, NZ 
CT=CO(I) *RW(I) 

JWT=QWZ*CO(0)  -D(0)/(Z(1)-Z(0))*(CO{1)-CO(0)  ) 

WRITEdO,  5000)  T,  Z(I),  CO  (I)  *1000000  .  dO ,  CT*1000000  .  dO , 
@JSF*10000.d0,  JWT*10000.d0 
PRINT*,  T,  Z{I),  CO(I) *1000000. do 
220  CONTINUE 

GOTO  400 
C 
C 

C  Boundary  Conditions 

************************************************************************ 

C 

C  LOADING  Phase;  C (Z=0, T=0 . .TFl) =0;  C (Z=L, T=0 . . TFl) =CI 
C  Note:  C  must  first  be  converted  from  ppm  (g/m'^3)  to  g/cm'^S) 

C 

Cl  300  CO(0) =CI/1000000,d0 

Cl  CO(NZ)=0.d0 

C2  300  CO(NZ)=0.d0 

Cl  C(0)=CI/1000000.d0 

C(NZ)=0.d0 
GOTO  500 
C 

C  UNLOADING  Phase:  C {Z=0 , T=TF1 . .TF2) =0;  C (Z=L, T=TF1 . . TF2) =0 
C 

400  CO(NZ)=0.d0 

Cl  C{0)=0.d0 
C{NZ) =0.d0 
GOTO  500 
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c 

c 

C  Solve  using  Thomas  algorithm 

************************************************************************ 

C 

500  G(l)  =E(1) /B(l) 

DO  520  J=2,NZ-1 

P(J)  =B(J)  -A(J)*G(J-1) 

G(J)  =E(J)/P(J) 

520  CONTINUE 

DO  600  N=1,NT 

C2  IF  (U.LE.NTL.AND.T+DELT.LE.200)  C (0) = (T+DELT) * (Cl 

C2  ©/lOOOOOO.dO) /200.d0 

C2  IF  (U.LE.NTL. AND. T+DELT. GT. 200)  C (0) =CI/1000000 . dO 
C2  IF  (U. GT.NTL. AND. T+DELT. LE. 22 00)  C (0) =CI/1000000 .dO 
C2  @* (1- (T-2000+DELT) /200.d0) 

C2  IF  (U. GT.NTL. AND. T+DELT. GT. 22 00)  C{0)=0.d0 

RHS  (1)  =-A(l)  *CO(0)  +  {2-B(l)  )*CO{l)  -E(l)  *CO{2)  -A(l)  *C(0) 
R(1)=RHS(1)/B(1) 

DO  540  J=2,NZ-2 

RHS ( J) =-A(J) *CO ( J-1) + (2-B (J) ) *CO (J) -E { J) *CO ( J+1) 

R ( J) = (RHS ( J) -A ( J) *R ( J-1) ) /P ( J) 

540  CONTINUE 

RHS (NZ-1) =-A(NZ-l) *CO (NZ-2) + (2-B (NZ-1) ) *CO (NZ-1) -E (NZ-1) 

@*CO(NZ) -E(NZ-1)*C(NZ) 

C (NZ-1) = (RHS (NZ-1) -A (NZ-1) *R(NZ-2) ) /P (NZ-1) 

DO  560  K=l,NZ-2 

C (NZ-l-K) =R (NZ-l-K) -G (NZ-l-K) *C (NZ-K) 

560  CONTINUE 

DO  580  1=0, NZ 
T=N*DELT+TO 
CT=C(I) *RW(I) 

JSF=QWZ*C (NZ) -D (NZ) / (Z (NZ) -Z (NZ-1) ) * (C (NZ) -C (NZ-1) ) 

JWT=QWZ*C(0) -D(0) / (Z(l) -Z(0) ) * (C(l) -C(0) ) 

C 

C  Convert  C  baclc  to  ppm  and  J  to  g/m^2  for  printing 
C 

WRITEdO,  5000)  T,  Z(I),  C  (I)  *1000000  .  dO ,  CT*1000000  .  dO , 
®JSF*10000.d0,  JWT*10000.d0 

PRINT*,  T,  Z(I),  C(I)*1000000,d0 
CO(I)=C(I) 

580  CONTINUE 

600  CONTINUE 

TO=T 

900  CONTINUE 
C 
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